library(clusterProfiler)
simpl_enrich <- simplify
library(enrichplot)
library(tidyverse)
library(ggpubr)
library(patchwork)
library(ggrepel)
source('00_util_scripts/mod_bplot.R')

proj.nm <- 'mission/FPP/zww_sa_mice/'

kegg_mva <-
  c('Acat1','Acat2','Hmgcs1','Hmgcs2','Hmgcr','Mvk','Pmvk','Mvd','Idi1','Idi2','Fdps')

mva_fct <- tibble(gene = kegg_mva, ordered = fct_inorder(kegg_mva))

kegg_d.mva <-
  c('FDFT1','LSS','SQLE','GGPS1','DOLK','PDSS1','PDSS2') |>
  str_to_title()

key_cytokine <- c('Il6','Ccl20','Tslp','Flt3l','Csf2','Tnf')

meta_kera <- read_csv('mission/FPP/zww_sa_mice/results/aureus_kera_meta.csv')

kera_infct_deg <-
  read_csv('mission/FPP/zww_sa_mice/results/kera.infect.deg.csv')

kera_v3_deg <-
  read_csv('mission/FPP/zww_sa_mice/results/mice_v3h_SA-PBS_deg.csv')

kera_v3l_deg <-
  read_csv('mission/FPP/zww_sa_mice/results/mice_v3l_SA-PBS_deg.csv')

sa_v3hvl_deg <-
  read_csv('mission/FPP/zww_sa_mice/results/mice_SA_v3h-v3l_deg.csv')

ct_v3hvl_deg <-
  read_csv('mission/FPP/zww_sa_mice/results/mice_ctrl_v3h-v3l_deg.csv')

# infect v3h vs ctrl v3h ------------
### FIG: SA v3h vs PBS v3h volcano ===========
kera_v3_deg |>
  filter(!str_detect(gene, 'Rik$|^Gm|AC116')) |>
  plot_bill_volc('SA infection','PBS control',
                 highlights = c(kegg_mva, key_cytokine)) +
  ggtitle('Infected Trpv3-hi-KC vs control Trpv3-hi-KC')

kera_v3_deg |>
  plot_pub_volc('SA Trpv3-hi KC','PBS Trpv3-hi KC', high_color = 'darkgreen',
                auto_n = 2)

g1 <- last_plot()

v3h_savpbs_goi <- kera_v3_deg |>
  filter(gene %in% c(kegg_mva, key_cytokine), avg_log2FC > 1, p_val_adj < .05)

g1 +
  geom_point(data = v3h_savpbs_goi, size = 1, color = 'purple') +
  geom_text_repel(data = v3h_savpbs_goi, aes(label = gene), color = 'black',
                  size = 2, box.padding = .1)

publish_source_plot('mice_v3h_SA-pbs_volc')

v3_inf_up <- kera_v3_deg |>
  filter(p_val_adj < .05 & avg_log2FC > 0)

v3_inf_down <- kera_v3_deg |>
  filter(p_val_adj < .05 & avg_log2FC < -2)

v3_inf_up_ora_gosim <- v3_inf_up$gene |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

v3_inf_up_ora_gosim |>
  pluck('result') |>
  head(12) |>
  plot_enrichment() +
  ggtitle('Upregulated GO BP pathway in Trpv3+ keratinocytes after infection')

v3_inf_down_ora <- v3_inf_down$gene |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

v3_inf_down_ora |>
  pluck('result') |>
  head(12) |>
  plot_enrichment('blue') +
  ggtitle('Downregulated GO BP pathway in Trpv3+ keratinocytes after infection')

v3_inf_up_ora_gosim@result |>
  as_tibble() |>
  filter(qvalue < .05) |>
  write_csv('mission/FPP/v3_inf-vs-pbs_up_go.csv')

v3_inf_down_ora@result |>
  as_tibble() |>
  filter(qvalue < .05) |>
  write_csv('mission/FPP/v3_inf-vs-pbs_down_go.csv')

### FIG: SA vs ctrl V3h-KC upgo ===========
v3_inf_up_ora_godf <-
  read_csv('mission/FPP/zww_sa_mice/results/v3_inf-vs-pbs_up_go.csv')

v3_inf_up_ora_godf |>
  publish_enrichment() +
  ggtitle('Upregulated GO BP pathway in V3-hi-KC: SA vs PBS')

publish_pdf('mission/FPP/micefig2/mice_v3h_SA-pbs_upgo.pdf', width = 60)

v3_inf_down_ora_godf <-
  read_csv('mission/FPP/v3_inf-vs-pbs_down_go.csv')

v3_inf_down_ora_godf |>
  head(12) |>
  plot_enrichment('blue') +
  ggtitle('Downregulated GO BP pathway in Trpv3+ keratinocytes after infection')

ggsave('mice_v3h_pbs-SA_downgo.pdf',
       width = 7, height = 4)

## interested genes in pathways? ---------
v3_inf_up_ora_godf |>
  separate_longer_delim(geneID, delim = '/') |>
  filter(geneID %in% chole_p_str) |>
  select(Description, qvalue)

v3_inf_up_ora_godf |>
  separate_longer_delim(geneID, delim = '/') |>
  filter(geneID %in% key_cytokine) |>
  select(Description, qvalue, geneID)

leuk_migr <- v3_inf_up_ora_godf |>
  separate_longer_delim(geneID, delim = '/') |>
  filter(Description == 'leukocyte migration') |>
  pull(geneID)

v3_inf_up |>
  filter(gene %in% leuk_migr)

sobj_kera |>
  filter(trpv3_status == 'Trpv3-high') |>
  ScaleData(features = leuk_migr) |>
  DoHeatmap(features = leuk_migr,
            group.by = 'orig.ident',
            angle = 0,
            vjust = 1,
            hjust = .5) +
  scale_fill_gradient2(high = 'red', low = 'blue')

g1 <- last_plot()

g1[[1]] + labs(title = 'GO BP pathway: Leukocyte migration')

ggsave('mice_v3h_SAup_GO1_heatmap.pdf',
       width = 6, height = 4)

small_catabo <- v3_inf_up_ora_godf |>
  separate_longer_delim(geneID, delim = '/') |>
  filter(Description == 'small molecule catabolic process') |>
  pull(geneID)

sobj_kera |>
  filter(trpv3_status == 'Trpv3-high') |>
  ScaleData(features = small_catabo) |>
  DoHeatmap(features = small_catabo,
            group.by = 'orig.ident',
            angle = 0,
            vjust = 1,
            hjust = .5) +
  scale_fill_gradient2(high = 'red', low = 'blue')

g2 <- last_plot()

g2[[1]] + labs(title = 'GO BP pathway: Small molecule catabolic process')

ggsave('mice_v3h_SAup_GO2_heatmap.pdf',
       width = 6, height = 4)

v3h.up.ros <- v3_inf_up_ora_godf |>
  filter(str_detect(Description, 'xygen')) |>
  separate_longer_delim(geneID, delim = '/') |>
  distinct(geneID) |>
  mutate(hit = anno_pmc_hits(str_glue('{geneID} ROS')))

v3h.up.ros

v3h.up.ros.hit <- v3h.up.ros |>
  arrange(desc(hit)) |>
  pull(geneID)

### GSEA GO ----------
v3h.savpbs.gslist <- kera_v3_deg |>
  filter(p_val_adj < .05) |>
  pull(avg_log2FC, name = 'gene') |>
  sort(decreasing = T)

v3h.savpbs.gsgo <- v3h.savpbs.gslist |>
  gseGO(OrgDb = 'org.Mm.eg.db',
        keyType = 'SYMBOL',
        ont = 'ALL', pvalueCutoff = 1)

v3h.savpbs.gsgo@result |>
  filter(NES > 0) |>
  plot_enrichment(metric = NES)

v3h.savpbs.gsgo@result |>
  as_tibble() |>
  filter(str_detect(Description, 'oxygen|stress|ndoplas')) |>
  select(ID,Description,NES,qvalue) |>
  DT::datatable()

v3h.savpbs.gsgo |>
  gseaplot2(
    'GO:2000379',
    title = 'GO pathway: positive regulation of ROS metabolic process\nV3h: SA vs PBS')

## ROS, ER stress hallmark genes --------
hllmk.ers <- c('HSPA5','DDIT3','EIF2A','EIF2AK3','ATF4')
kera_v3l_deg |>
  filter(gene %in% str_to_title(hllmk.ers))

v3h.upers <- v3_inf_up_ora_gosim@result |>
  as_tibble() |>
  filter(str_detect(Description, 'reticulum stress')) |> select(1,2,Count)
  separate_longer_delim(geneID, '/') |>
  distinct(geneID) 

v3h.upers.hit <- v3h.upers |>
  mutate(hit = anno_pmc_hits(str_glue('ER stress {geneID}')))
  
v3h.upers.hit <- v3h.upers.hit |>
  mutate(gene = geneID, .keep = 'unused') |>
  left_join(v3_inf_up)

v3h.upers.hit |>
  write_csv('mission/FPP/zww_sa_mice/results/v3h.savpbs.up.ERs.hit.csv')

v3h.upers.hit |>
  ggplot(aes(pct.1, hit, label = gene)) +
  geom_point()

plotly::ggplotly()

v3h.upers.hit9 <-
  tibble(gene = v3h.upers, hit = v3h.upers.hit) |>
  left_join(v3_inf_up) |>
  slice_max(hit, n = 9) |>
  pull(gene)

v3h.upers.top9 <- v3_inf_up |>
  filter(gene %in% v3h.upers, pct.1 > .39) |>
  slice_max(avg_log2FC, n = 9) |>
  pull(gene)

kera_v3l_deg |>
  filter(gene %in% v3h.upers.top9)

v3h.upapop <-
  v3_inf_up_ora_gosim@result |>
  as_tibble() |>
  filter(str_detect(Description, 'intrinsic apop')) |>
  separate_longer_delim(geneID, '/') |>
  distinct(geneID) |>
  pull(geneID)

v3h.upapop.hit <- v3h.upapop |>
  map(str_c, ' apoptosis') |>
  anno_pmc_hits()

v3h.upapop.hit30 <-
  tibble(gene = v3h.upapop, hit = v3h.upapop.hit) |>
  left_join(v3_inf_up) |>
  slice_max(hit, n = 30) |>
  pull(gene)

### custom gsea all cell types -------
savpbs.skin.logfc <-
  read_csv('mission/FPP/zww_sa_mice/results/savpbs.skin.all.logfc.csv')

savpbs.skin.logfc |>
  filter(gene == 'Cxcl12')

#### UPR chaperone -------
upers.chpr <- read_csv('mission/FPP/UPR.associated.chaperone.gene.csv')

savpbs.skin.logfc |>
  filter(gene %in% str_to_title(upers.chpr$SYMBOL)) |>
  write_csv('fig.E3.ERS.chaperone.skin.SA.vs.PBS.v3h.v3l.csv')

sa.fine <- savpbs.skin.logfc$cluster |> unique()

skin.savpbs.ers <- sa.fine |>
  map(\(x){savpbs.skin.logfc |>
      filter(cluster == x) |>
      pull(avg_log2FC, name = gene) |>
      sort(decreasing = T) |>
      gseGO(ont = 'ALL', OrgDb = 'org.Mm.eg.db',
            keyType = 'SYMBOL',pvalueCutoff = 1) |>
      pluck('result')}) |>
  set_names(sa.fine) |>
  list_rbind(names_to = 'cell.type') |>
  as_tibble()

skin.savpbs.ers |>
  select(-c(setSize, enrichmentScore)) |>
  write_csv('mission/FPP/zww_sa_mice/results/skin.all.savpbs.go.gsea.csv')

skin.savpbs.ers <-
  read_csv('mission/FPP/zww_sa_mice/results/skin.all.savpbs.go.gsea.csv')

plot.gsea.dot <- function(df) {
  df |> ggplot(aes(.data$cell.type, str_wrap(.data$Description, 50),
             size = -log10(.data$pvalue), color = .data$NES)) +
    geom_point() +
    scale_size(range = c(0,3)) +
    scale_color_distiller(palette = 'RdYlBu') +
    theme_jpub +
    rotate_x_text(45) +
    labs(title = 'Differential pathway enrichment in SA vs PBS',
         y = 'GO term', x = 'Cell type')
}

skin.savpbs.ers |>
  filter(str_detect(Description, 'reactive oxygen')) |>
  plot.gsea.dot()

publish_pdf('mission/FPP/micefig2/skin.savpbs.ros.gsea.pdf', width = 70)

skin.savpbs.ers |>
  filter(str_detect(Description, 'TOR')) |>
  write_csv('fig.H.TOR.signaling.in.SA.vs.PBS.skin.csv')

skin.savpbs.ers |>
  filter(str_detect(Description, 'TOR')) |>
  plot.gsea.dot()

publish_pdf('mission/FPP/micefig2/skin.savpbs.tor.gsea.pdf', width = 60)

skin.savpbs.ers |>
  filter(str_detect(Description, '[^o]autophag')) |>
  plot.gsea.dot()

publish_pdf('mission/FPP/micefig2/skin.savpbs.autophagy.gsea.pdf', width = 70)

skin.savpbs.ers |>
  filter(str_detect(Description, 'pyropto|inflammasom')) |>
  plot.gsea.dot()

publish_pdf('mission/FPP/micefig2/skin.savpbs.inflammasome.gsea.pdf', width = 75)

skin.savpbs.ers |>
  filter(str_detect(Description, 'toll')) |>
  plot.gsea.dot()

publish_pdf('mission/FPP/micefig2/skin.savpbs.TLR.gsea.pdf', width = 70)

skin.savpbs.ers |>
  filter(str_detect(Description, 'interleukin-(1$|1 )')) |>
  plot.gsea.dot()

publish_pdf('mission/FPP/micefig2/skin.savpbs.IL1.gsea.pdf', width = 70)

skin.savpbs.ers |>
  filter(str_detect(Description, 'epidermal growth factor')) |>
  write_csv('fig.F.EGF.signaling.in.SA.vs.PBS.skin.csv')

skin.savpbs.ers |>
  filter(str_detect(Description, 'epidermal growth factor')) |>
  plot.gsea.dot()

publish_pdf('mission/FPP/micefig2/skin.savpbs.EGF.gsea.pdf', width = 70)

### Overview GO dotplot ----------
#### SA vs PBS ---------
marker.go <- skin.savpbs.ers |>
  filter(ONTOLOGY == 'BP') |>
  slice_max(NES, n = 2, by = cell.type, with_ties = F)
  
ordered.go <- marker.go |>  
  arrange(cell.type) |>
  mutate(term = fct_inorder(Description), .keep = 'used')

skin.savpbs.ers |>
  filter(Description %in% marker.go$Description) |>
  left_join(ordered.go) |>
  ggplot(aes(.data$cell.type, term,
             size = -log10(.data$pvalue), color = .data$NES)) +
  geom_point() +
  scale_size(range = c(0,3)) +
  scale_color_distiller(palette = 'RdYlBu') +
  theme_jpub +
  rotate_x_text(45) +
  scale_y_discrete(label = \(x)str_wrap(x, 45)) +
  labs(title = 'Differential pathway enrichment in SA vs PBS',
       y = 'GO term', x = 'Cell type')

publish_pdf('mission/FPP/micefig2/skin.savpbs.allrep.gsea.pdf', width = 80,
            height = 80)

last_plot() |>
  pluck('data') |>
  write_csv('mission/FPP/pub_source_data/R3.Q10.figC.SA.vs.PBS.skin.GO.csv')

##### T & T subsets ----------
skin.savpbs.ers |>
  filter(cell.type == 'T cells') |>
  mutate(direction = ifelse(NES > 0, 'Up', 'Down')) |>
  slice_min(p.adjust, n = 5, with_ties = F, by = direction) |>
  mutate(Description = str_wrap(Description, 40) |>
           fct_reorder(NES)) |>
  ggplot(aes(Description, NES, fill = p.adjust)) +
  geom_col() +
  coord_flip() +
  theme_pubr(legend = 'right') +
  labs_pubr() +
  scale_fill_gradient(low = 'red', high = 'black') +
  labs(title = 'GO pathway enrichment of T cells: SA vs PBS')


#### Within SA ----------
sa.skin.allmark <-
  read_csv('mission/FPP/zww_sa_mice/results/sa.skin.all.marker.csv')

sa.skin.allgo <- sa.fine |>
  map(\(x){sa.skin.allmark |>
      filter(cluster == x) |>
      pull(avg_log2FC, name = gene) |>
      sort(decreasing = T) |>
      gseGO(ont = 'ALL', OrgDb = 'org.Mm.eg.db',
            keyType = 'SYMBOL',pvalueCutoff = 1) |>
      pluck('result')}) |>
  set_names(sa.fine) |>
  list_rbind(names_to = 'cell.type') |>
  as_tibble()

sa.skin.allgo |>
  select(-c(setSize, enrichmentScore)) |>
  write_csv('mission/FPP/zww_sa_mice/results/skin.all.sa.go.gsea.csv')

sa.skin.allgo <-
  read_csv('mission/FPP/zww_sa_mice/results/skin.all.sa.go.gsea.csv')

sa.go.marker <- sa.skin.allgo |>
  filter(ONTOLOGY == 'BP') |>
  slice_max(NES, n = 2, by = cell.type, with_ties = F)

sa.go.ordered <- sa.go.marker |>  
  arrange(cell.type) |>
  mutate(term = fct_inorder(Description), .keep = 'used')

sa.skin.allgo |>
  filter(Description %in% sa.go.marker$Description) |>
  left_join(sa.go.ordered) |>
  ggplot(aes(.data$cell.type, term,
             size = -log10(.data$pvalue), color = .data$NES)) +
  geom_point() +
  scale_size(range = c(0,3)) +
  scale_color_distiller(palette = 'RdYlBu') +
  theme_jpub +
  rotate_x_text(45) +
  scale_y_discrete(label = \(x)str_wrap(x, 45)) +
  labs(title = 'Differential pathway enrichment in SA skin',
       y = 'GO term', x = 'Cell type')

publish_pdf('mission/FPP/micefig2/skin.sa.allrep.gsea.pdf', width = 80,
            height = 80)

#### Within PBS ----------
pbs.skin.allmark <-
  read_csv('mission/FPP/zww_sa_mice/results/pbs.skin.all.marker.csv')

pbs.skin.allgo <- sa.fine |>
  map(\(x){pbs.skin.allmark |>
      filter(cluster == x) |>
      pull(avg_log2FC, name = gene) |>
      sort(decreasing = T) |>
      gseGO(ont = 'ALL', OrgDb = 'org.Mm.eg.db',
            keyType = 'SYMBOL',pvalueCutoff = 1) |>
      pluck('result')}) |>
  set_names(sa.fine) |>
  list_rbind(names_to = 'cell.type') |>
  as_tibble()

pbs.skin.allgo |>
  select(-c(setSize, enrichmentScore)) |>
  write_csv('mission/FPP/zww_sa_mice/results/skin.all.pbs.go.gsea.csv')

pbs.skin.allgo <-
  read_csv('mission/FPP/zww_sa_mice/results/skin.all.pbs.go.gsea.csv')

pbs.go.marker <- pbs.skin.allgo |>
  filter(ONTOLOGY == 'BP', str_detect(Description, 'synapse', negate = T)) |>
  slice_max(NES, n = 2, by = cell.type, with_ties = F)

pbs.go.ordered <- pbs.go.marker |>  
  arrange(cell.type) |>
  mutate(term = fct_inorder(Description), .keep = 'used')

pbs.skin.allgo |>
  filter(Description %in% pbs.go.marker$Description) |>
  left_join(pbs.go.ordered) |>
  ggplot(aes(.data$cell.type, term,
             size = -log10(.data$pvalue), color = .data$NES)) +
  geom_point() +
  scale_size(range = c(0,3)) +
  scale_color_distiller(palette = 'RdYlBu') +
  theme_jpub +
  rotate_x_text(45) +
  scale_y_discrete(label = \(x)str_wrap(x, 40)) +
  labs(title = 'Differential pathway enrichment in PBS skin',
       y = 'GO term', x = 'Cell type')

publish_pdf('mission/FPP/micefig2/skin.pbs.allrep.gsea.pdf', width = 80,
            height = 80)

### Arachidonic acid synthesis -------
arachi.syn <-
'18778(Pla2g1b) 18780(Pla2g2a) 18781(Pla2g2c) 18782(Pla2g2d) 18784(Pla2g5) 237625(Pla2g3) 26565(Pla2g10) 26970(Pla2g2e) 26971(Pla2g2f) 66350(Pla2g12a) 69836(Pla2g12b)' |>
  str_extract_all('Pla\\w+') |>
  pluck(1)

kera_v3_deg |>
  filter(str_starts(gene, 'Pla2g4|Plcb1'))

kera_v3l_deg |>
  filter(gene %in% arachi.syn)

# infect v3l vs ctrl v3l ------------
### FIG: SA v3l vs PBS v3l kera volcano ===========
kera_v3l_deg |>
  filter(!str_detect(gene, 'Rik$|^Gm|AC116'), p_val_adj != 0) |>
  plot_bill_volc('SA V3-lo KC','PBS V3-lo KC',
                 highlights = c(kegg_mva, key_cytokine))

kera_v3l_deg |>
  filter(p_val_adj != 0) |>
  plot_pub_volc('SA V3-lo KC','PBS V3-lo KC',
                 highlights = c(kegg_mva, key_cytokine))

publish_pdf('mission/FPP/micefig2/mice_v3l_SA-pbs_volc.pdf')

v3l_inf_up <- kera_v3l_deg |>
  filter(p_val_adj < .05 & avg_log2FC > 1)

v3l_inf_down <- kera_v3l_deg |>
  filter(p_val_adj < .05 & avg_log2FC < -1)

v3l_inf_upgo <- v3l_inf_up$gene |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

v3l_inf_upgo@result |>
  plot_enrichment() +
  ggtitle('Upregulated GO BP pathway in Trpv3-lo keratinocytes after infection')

v3l_inf_upgo@result |>
  publish_enrichment() +
  ggtitle('Upregulated GO BP pathway in V3-lo KC: SA vs PBS')

publish_pdf('mission/FPP/micefig2/v3l_sa_vs_pbs_upgo.pdf', width = 65)

v3_inf_down_ora <- v3_inf_down$gene |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

v3_inf_down_ora |>
  pluck('result') |>
  head(12) |>
  plot_enrichment('blue') +
  ggtitle('Downregulated GO BP pathway in Trpv3+ keratinocytes after infection')

v3_inf_up_ora_gosim@result |>
  as_tibble() |>
  filter(qvalue < .05) |>
  write_csv('mission/FPP/v3_inf-vs-pbs_up_go.csv')

v3_inf_down_ora@result |>
  as_tibble() |>
  filter(qvalue < .05) |>
  write_csv('mission/FPP/v3_inf-vs-pbs_down_go.csv')

### FIG: SA vs ctrl V3h-KC upgo ===========
v3_inf_up_ora_godf <-
  read_csv('mission/FPP/zww_sa_mice/results/v3_inf-vs-pbs_up_go.csv')

v3_inf_up_ora_godf |>
  head(10) |>
  plot_enrichment() +
  ggtitle('Upregulated GO BP pathway in SA-infected vs control V3-hi-KC')

publish_pdf('micefig/mice_v3h_SA-pbs_upgo.pdf', width = 60)

v3_inf_down_ora_godf <-
  read_csv('mission/FPP/v3_inf-vs-pbs_down_go.csv')

v3_inf_down_ora_godf |>
  head(12) |>
  plot_enrichment('blue') +
  ggtitle('Downregulated GO BP pathway in Trpv3+ keratinocytes after infection')

### FIG: v3l SA-PBS upgo ===============
v3l_sa_upgo <- read_csv('mission/FPP/v3l_kera_SA-PBS_deg.csv') |>
  filter(p_val_adj < .05 & avg_log2FC > 2) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

v3l_sa_upgo@result |>
  head(10) |>
  publish_enrichment()

publish_pdf('micefig/mice_SA_v3l_SA-PBS_upgo.pdf', width = 60)

read_csv('mission/FPP/inf_v3h-vs-v3l_down_go.csv') |>
  head(12) |>
  plot_enrichment('blue') +
  ggtitle('Downregulated GO BP pathway in Trpv3+ vs Trpv3- keratinocytes')

ggsave('mice_SA_kera_v3h-v3l_downgo.pdf',
       width = 7, height = 4)

v3_mrkr_up_go@result |>
  as_tibble() |>
  filter(str_detect(Description, 'retinoid meta')) |>
  pull(geneID)

### GSEA GO ----------
v3l.savpbs.gslist <- kera_v3l_deg |>
  filter(p_val_adj < .05) |>
  pull(avg_log2FC, name = 'gene') |>
  sort(decreasing = T)

v3l.savpbs.gsgo <- v3l.savpbs.gslist |>
  gseGO(OrgDb = 'org.Mm.eg.db',
        keyType = 'SYMBOL',
        ont = 'ALL', pvalueCutoff = 1)

v3l.savpbs.gsgo@result |>
  as_tibble() |>
  filter(str_detect(Description, 'oxygen|stress|ndoplas'), qvalue < .1) |>
  DT::datatable()

v3l.savpbs.gsgo |>
  gseaplot2(
    'GO:2000379',
    title = 'GO pathway: positive regulation of ROS metabolic process\nV3l: SA vs PBS')

# infect kera vs ctrl kera -------
### FIG: SA vs ctrl kera volcano ===========
kera_infct_deg |>
  filter(gene %in% c(key_cytokine, kegg_mva))

kera_infct_deg |>
  filter(!str_detect(gene, 'Rik$|^Gm')) |>
  plot_pub_volc('keratinocytes', highlights = c(key_cytokine, kegg_mva))

publish_pdf('micefig/mice_pbs-kera_SA-kera_volcano.pdf')

### FIG: SA vs ctrl kera upgo ===========
kr_inf_up <- kera_infct_deg |>
  filter(p_val_adj < .05 & avg_log2FC > 2)

kr_inf_down <- kera_infct_deg |>
  filter(p_val_adj < .05 & avg_log2FC < -2)

kr_inf_up_ora_go <- kr_inf_up$gene |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

kr_inf_up_ora_go |>
  pluck('result') |>
  head(12) |>
  plot_enrichment() +
  ggtitle('Upregulated GO BP pathway in all keratinocytes after infection')

kr_inf_down_ora <- kr_inf_down$gene |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

kr_inf_down_ora |>
  pluck('result') |>
  head(12) |>
  plot_enrichment('blue') +
  ggtitle('Downregulated GO BP pathway in all keratinocytes after infection')

kr_inf_up_ora_go@result |>
  as_tibble() |>
  filter(qvalue < .05) |>
  write_csv('mission/FPP/kr_inf-vs-pbs_up_go.csv')

kr_inf_down_ora@result |>
  as_tibble() |>
  filter(qvalue < .05) |>
  write_csv('mission/FPP/kr_inf-vs-pbs_down_go.csv')

read_csv('mission/FPP/kr_inf-vs-pbs_up_go.csv') |>
  head(10) |>
  publish_enrichment()

publish_pdf('micefig/mice_all_kera_pbs-SA_upgo.pdf', width = 65)

read_csv('mission/FPP/kr_inf-vs-pbs_down_go.csv') |>
  head(12) |>
  plot_enrichment('blue') +
  ggtitle('Downregulated GO BP pathway in all keratinocytes after infection')

ggsave('mice_all_kera_pbs-SA_downgo.pdf',
       width = 7, height = 4)

### GSEA GO ----------
kc.savpbs.gslist <- kera_infct_deg |>
  filter(p_val_adj < .05) |>
  pull(avg_log2FC, name = 'gene') |>
  sort(decreasing = T)

kc.savpbs.gsgo <- kc.savpbs.gslist |>
  gseGO(OrgDb = 'org.Mm.eg.db',
        keyType = 'SYMBOL',
        ont = 'ALL', pvalueCutoff = 1)

kc.savpbs.gsgo@result |>
  as_tibble() |>
  filter(str_detect(Description, 'oxygen|stress|ndoplas'), qvalue < .1) |>
  DT::datatable()

kc.savpbs.gsgo |>
  gseaplot2(
    'GO:2000379',
    title = 'GO pathway: positive regulation of ROS metabolic process\nKC: SA vs PBS')

sa_kc_upr_nes <- kc.savpbs.gsgo@result |>
  as_tibble() |>
  filter(str_detect(Description, 'unfold')) |>
  select(ID, Description, NES, qvalue, leading_edge)

sa_kc_upr_nes

# SA: v3h vs v3l ---------
## FIG: SA v3h vs v3l volcano =============
sa_v3hvl_deg |>
  filter(!str_detect(gene, 'AC116') & p_val_adj != 0) |>
  plot_bill_volc('V3-hi KC',group2 = 'V3-lo KC',
                 highlights = c(key_cytokine, 'Trpv3'))

sa_v3hvl_deg |>
  filter(!str_detect(gene, 'AC116') & p_val_adj != 0) |>
  plot_pub_volc('SA Trpv3-hi KC', group2 = 'SA Trpv3-lo KC',
                high_color = 'darkgreen', auto_n = 2)

g1 <- last_plot()

sa_v3hvl_goi <- sa_v3hvl_deg |>
  filter(gene %in% c(key_cytokine, 'Mvk', 'Pmvk', 'Trpv3'),
         avg_log2FC > 1, p_val_adj < .05)

g1 +
  geom_point(data = sa_v3hvl_goi, size = 1, color = 'purple') +
  geom_text_repel(data = sa_v3hvl_goi, aes(label = gene), color = 'black',
                  size = 2, box.padding = .1)

publish_pdf('mice_SA_v3h-v3l_volcano.pdf')

inf_v3hvl_cyto <- trpv3_kera_marker |>
  filter(gene %in% key_cytokine)

inf_v3hvl_cyto |>
  ggplot(aes(x = 'Trpv3-hi vs Trpv3-lo', y = gene, color = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point() +
  scale_color_gradient2(high = 'red',low = 'blue') +
  theme_pubr(legend = 'right') +
  labs_pubr() +
  labs(x = '')

ggsave('mice_SA_v3h-v3l-logfc_cytokine_bubbleplot.pdf',
       width = 4, height = 3.5)

## GO ORA ---------
v3_mrkr_up <- sa_v3hvl_deg |>
  filter(p_val_adj < .05 & avg_log2FC > 1)

v3_mrkr_down <- sa_v3hvl_deg |>
  filter(p_val_adj < .05 & avg_log2FC < -1)

v3_mrkr_up_go <- v3_mrkr_up$gene |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

v3_mrkr_down_go <- v3_mrkr_down$gene |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

v3_mrkr_up_go@result |>
  publish_enrichment() +
  ggtitle('Upregulated GO BP pathway in SA group: V3-hi KC vs V3-lo KC')

publish_pdf('mission/FPP/micefig2/mice_SA_v3hvl_upgo.pdf', width = 65)

v3_mrkr_down_go@result |>
  head(12) |>
  plot_enrichment('blue') +
  ggtitle('Downregulated GO BP pathway in Trpv3+ vs Trpv3- keratinocytes')

v3_mrkr_up_go@result |>
  as_tibble() |>
  filter(qvalue < .05) |>
  write_csv('mission/FPP/inf_v3h-vs-v3l_up_go.csv')

v3_mrkr_down_go@result |>
  as_tibble() |>
  filter(qvalue < .05) |>
  write_csv('mission/FPP/inf_v3h-vs-v3l_down_go.csv')

### FIG: SA v3h vs v3l upgo ============
read_csv('mission/FPP/inf_v3h-vs-v3l_up_go.csv') |>
  head(10) |>
  publish_enrichment()

publish_pdf('micefig/mice_SA_kera_v3h-v3l_upgo.pdf', width = 60)

### GSEA GO ----------
sa.v3hvl.gslist <- sa_v3hvl_deg |>
  filter(p_val_adj < .05) |>
  pull(avg_log2FC, name = 'gene') |>
  sort(decreasing = T)

sa.v3hvl.gsgo <- sa.v3hvl.gslist |>
  gseGO(OrgDb = 'org.Mm.eg.db',
        keyType = 'SYMBOL',
        ont = 'ALL', pvalueCutoff = 1, eps = 0, nPermSimple = 10000)

sa.v3hvl.gsgo@result |>
  as_tibble() |>
  filter(str_detect(Description, 'apopto'), qvalue < .05) |>
  DT::datatable()

nperm1000 <- sa.v3hvl.gsgo@result$ID

sa.v3hvl.gsgo <- sa.v3hvl.gsgo |>
  simplify()

sa.v3hvl.gsgo@result |>
  as_tibble() |>
  filter(!(ID %in% nperm1000)) |>
  DT::datatable()

sa.v3hvl.gsgo |>
  gseaplot2(
    'GO:0072593',
    title = 'GO pathway: ROS metabolic process\nSA: V3h vs V3l')

sa.v3hvl.gsgo |>
  gseaplot2(
    'GO:0097191',
    title = 'GO pathway: extrinsic apoptotic signaling pathway\nSA: V3h vs V3l')

## Bcl2/Bax apoptosis ------------
bb.apop <-
  c('Bcl2','Bax','Bak1','Diablo','Htra2','Cyc','Apaf1','Fadd','Casp3','Tnfrsf1a')

sa_v3hvl_deg |>
  filter(gene %in% bb.apop)

search_go_term('apoptosis')

# PBS: v3h vs v3l --------------
## FIG: pbs v3h vs v3l volcano =============
ct_v3hvl_deg |>
  filter(!str_detect(gene, 'AC116') & p_val_adj != 0) |>
  plot_bill_volc('PBS Trpv3-hi KC',group2 = 'PBS Trpv3-lo KC',
                 highlights = c(key_cytokine,kegg_mva)) +
  ggtitle('PBS Trpv3-hi KC vs PBS Trpv3-lo KC')

ct_v3hvl_deg |>
  filter(!str_detect(gene, 'AC116') & p_val_adj != 0) |>
  plot_pub_volc('PBS V3-hi KC',group2 = 'PBS V3-lo KC',
                 highlights = c(key_cytokine,kegg_mva))

publish_pdf('mission/FPP/micefig2/mice_pbs_v3h-v3l_volcano.pdf')

## GO ORA ---------
ct_v3hvl_up <- ct_v3hvl_deg |>
  filter(p_val_adj < .05 & avg_log2FC > 1)

ct_v3hvl_down <- ct_v3hvl_deg |>
  filter(p_val_adj < .05 & avg_log2FC < -1)

ct_v3hvl_upgo <- ct_v3hvl_up$gene |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

v3_mrkr_down_go <- v3_mrkr_down$gene |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

ct_v3hvl_upgo@result |>
  plot_enrichment() +
  ggtitle('Upregulated GO BP pathway in control group Trpv3-hi vs Trpv3-lo KC')

ct_v3hvl_upgo@result |>
  publish_enrichment() +
  ggtitle('Upregulated GO BP pathway in PBS group: V3-hi KC vs V3-lo KC')

publish_pdf('mission/FPP/micefig2/pbs_v3hvl_upgo.pdf', width = 65)

ct_v3hvl_upgo |>
  as_tibble() |>
  filter(qvalue < .05) |>
  write_csv('mission/FPP/zww_sa_mice/pbs_v3hvl_upgo.csv')

read_tsv('mission/FPP/zww_sa_mice/results/pbs.v3hvl.upgo.digest.tsv') |>
  plot_enrichment() +
  ggtitle('Upregulated GO pathway in PBS group v3h vs v3l KC')

v3_mrkr_down_go@result |>
  as_tibble() |>
  filter(qvalue < .05) |>
  write_csv('mission/FPP/inf_v3h-vs-v3l_down_go.csv')

## GSEA GO ----------
pbs.v3hvl.gslist <- ct_v3hvl_deg |>
  filter(p_val_adj < .05) |>
  pull(avg_log2FC, name = 'gene') |>
  sort(decreasing = T)

pbs.v3hvl.gsgo <- pbs.v3hvl.gslist |>
  gseGO(OrgDb = 'org.Mm.eg.db',
        keyType = 'SYMBOL',
        ont = 'ALL', pvalueCutoff = 1)

pbs.v3hvl.gsgo@result |>
  as_tibble() |>
  filter(str_detect(Description, 'oxygen|stress|ndoplas'), qvalue < .1) |>
  DT::datatable()

pbs.v3hvl.gsgo |>
  gseaplot2(
    'GO:0072593',
    title = 'GO pathway: ROS metabolic process\nPBS: V3h vs V3l')

### GSEA KEGG --------
pbs.v3hvl.deg <-
ct_v3hvl_deg$gene |>
  bitr(fromType = 'SYMBOL', toType = 'ENTREZID', OrgDb = 'org.Mm.eg.db') |>
  left_join(ct_v3hvl_deg, join_by(SYMBOL == gene))

pbs.v3hvl.kelist <- pbs.v3hvl.deg |>
  filter(p_val_adj < .05) |>
  pull(avg_log2FC, name = ENTREZID) |>
  sort(decreasing = T)

pbs.v3hvl.gske <- pbs.v3hvl.kelist |>
  gseKEGG(organism = 'mmu')

pbs.v3hvl.gske@result |>
  as_tibble() |>
  filter(str_detect(Description, 'xygen|tress|ndoplas'), qvalue < .1) |>
  DT::datatable()

pbs.v3hvl.gsgo |>
  gseaplot2(
    'GO:0072593',
    title = 'GO pathway: ROS metabolic process\nPBS: V3h vs V3l')

pbs.v3hvl.gspc <- pbs.v3hvl.gslist |>
  gsePC()

pbs.v3hvl.gswp@result |>
  as_tibble() |>
  #filter(str_detect(Description, 'xygen|tress|ndoplas'), qvalue < .1) |>
  DT::datatable()

# compare cell frac change upon infection --------
stat_kera <- meta_kera |>
  count(orig.ident, seurat_clusters) |>
  group_by(orig.ident) |>
  calc_frac_conf_on_grouped_count()

pval_kera <- meta_kera |>
  test_on_grouped_count(orig.ident, seurat_clusters) |>
  ungroup() |>
  mutate(p.adj = p.adjust(p.value, method = 'BH'),
         seurat_clusters = subtype)

stat_kera |>
  ungroup() |>
  select(orig.ident, seurat_clusters, fraction) |>
  pivot_wider(names_from = orig.ident, values_from = fraction) |>
  left_join(pval_kera) |>
  mutate(fraction_log2fc = log2(infected / PBS),
         type = case_when(p.adj < .05 & fraction_log2fc < 0 ~ 'Decreased after infection',
                          p.adj < .05 & fraction_log2fc > 0 ~ 'Increased after infection',
                          .default = 'NS'),
         clusters = as.character(seurat_clusters) |> fct_reorder(seurat_clusters)) |>
  ggplot(aes(clusters, fraction_log2fc, fill = type)) +
  geom_col() +
  coord_flip() +
  theme_pubr(legend = 'right') +
  scale_fill_manual(values = c('blue','red','grey')) +
  labs(y = 'Log2 fold changes after infection',
       title = 'Cell fraction changes of keratinocytes after infection')

## FIG: v3h frac in PBS & SA ===========
v3_frac_conf <- meta_kera |>
  calc_frac_conf_on_grouped_count(orig.ident, trpv3_status)

v3_frac_pval <- meta_kera |>
  test_on_grouped_count(orig.ident, trpv3_status)

v3_frac_conf |>
  filter(trpv3_status == 'Trpv3-high') |>
  mutate(group = fct_relevel(orig.ident, 'PBS')) |>
  ggplot(aes(group, fraction, ymin = conf.low, ymax = conf.high, fill = group)) +
  geom_col() +
  geom_errorbar(width = .5) +
  scale_fill_manual(values = c('blue','red'), labels = c('PBS','SA')) +
  labs(y = 'Fraction in all keratinocytes', x = 'Group',
       title = 'Fraction of Trpv3-hi KC in mice skin') +
  theme_pubr(legend = 'right') +
  theme_jpub

publish_pdf('mission/FPP/micefig2/mice_v3h_krea_SA-PBS_fraction_barplot.pdf',
            width = 40)

v3_frac_conf |>
  filter(trpv3_status == 'Trpv3-low') |>
  mutate(group = fct_relevel(orig.ident, 'PBS')) |>
  ggplot(aes(group, fraction, ymin = conf.low, ymax = conf.high, fill = group)) +
  geom_col() +
  geom_errorbar(width = .5) +
  scale_fill_manual(values = c('blue','red'), labels = c('PBS','SA')) +
  labs(y = 'Fraction in all keratinocytes', x = 'Group',
       title = 'Fraction of Trpv3-low KC in mice skin') +
  theme_pubr(legend = 'right') +
  theme_jpub

publish_pdf('mission/FPP/micefig3/mice_v3l_krea_SA-PBS_fraction_barplot.pdf',
            width = 40)

v3_frac_conf |>
  write_csv('mission/FPP/pub_source_data/v3h.v3l.in.KC.fraction.csv')

## FIG: frac of other immune cells ---------
skin.meta <- read_csv('mission/FPP/zww_sa_mice/results/aureus_skin_meta.csv')

skin.meta |>
  filter(str_detect(manual_main, 'Mac|NK|T|Gra')) |>
  calc_frac_conf_on_grouped_count(orig.ident, manual_main) |>
  mutate(group = fct_relevel(orig.ident, 'PBS')) |>
  ggplot(aes(group, fraction,
             ymin = conf.low, ymax = conf.high, color = group)) +
  geom_col(fill = 'white') +
  geom_errorbar(width = .5) +
  facet_wrap(~manual_main, scales = 'free_y') +
  scale_color_hue(direction = -1) +
  labs(y = 'Fraction in all immune cells', x = 'Group', color = 'Group',
       title = 'Fraction of immune cells in mice skin') +
  theme_pubr(legend = 'right') +
  theme_jpub

publish_pdf('mission/FPP/micefig2/skin.immune.frac.pdf', width = 60)

# compare 8 vs 10 kera ----------
kera8v10 <- read_csv('mission/FPP/zww_sa_mice/kera8v10.deg.csv')

key_signaling <- c('Ptk2b','Kras','Nfkb1','Nfkb2','Calm1','Ppp3cc','Nfam1')

kera8v10 |>
  filter(gene %in% key_signaling)

## go ora ------------
kera8_upgo <- kera8v10 |>
  filter(avg_log2FC > 2 & p_val_adj < .05) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           qvalueCutoff = .05,
           readable = T)
  
kera8_upgo@result |>
  head(12) |>
  plot_enrichment() +
  ggtitle('Upregulated GO pathways in cluster 8 keratinocytes')

kera8_upgo@result |>
  as_tibble() |>
  filter(p.adjust < .05) |>
  DT::datatable()

kera10_upgo <- kera8v10 |>
  filter(avg_log2FC < 2 & p_val_adj < .05) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           qvalueCutoff = .05,
           readable = T)

kera10_upgo@result |>
  head(12) |>
  plot_enrichment() +
  ggtitle('Upregulated GO pathways in cluster 10 keratinocytes')

## go gsea ---------------
sig8v10 <- kera8v10 |>
  filter(p_val_adj < .05)

gsgo8v10 <- sig8v10$avg_log2FC |>
  set_names(sig8v10$gene) |>
  sort(decreasing = T) |>
  gseGO(OrgDb = 'org.Mm.eg.db',
        keyType = 'SYMBOL',
        ont = 'BP',
        pvalueCutoff = 1)

gsgo8v10@result |>
  as_tibble() |>
  filter(str_detect(Description, 'alcium')) |>
  DT::datatable()

gsgo8v10 |>
  gseaplot2(geneSetID = c('GO:0071277','GO:0019722'), pvalue_table = T)

gsgo10v8 <- -sig8v10$avg_log2FC |>
  set_names(sig8v10$gene) |>
  sort(decreasing = T) |>
  gseGO(OrgDb = 'org.Mm.eg.db',
        keyType = 'SYMBOL',
        ont = 'BP')

gsgo10v8@result |>
  as_tibble() |>
  slice_max(NES, n =3)

# pathview -----
v3h.sa.vs.pbs <- kera_v3_deg |>
  pull(avg_log2FC, name = gene)
  
v3l.sa.vs.pbs <- kera_v3l_deg |>
  pull(avg_log2FC, name = gene)

kera_v3_deg |>
  filter(gene %in% kegg_d.mva) |>
  mutate(gene = fct_relevel(gene, 'Fdft1', 'Sqle','Lss') |> fct_rev()) |>
  ggplot(aes('', gene, fill = avg_log2FC)) +
  geom_tile() +
  scale_fill_gradient2(low = 'blue', high = 'red', limits = c(-1,2),
                       na.value = 'grey80') +
  theme_pubr(legend = 'right') +
  labs(x = '', y = 'Gene')

kera_v3_deg |>
  right_join(mva_fct) |>
  ggplot(aes('', fct_rev(ordered), fill = avg_log2FC)) +
  geom_tile() +
  scale_fill_gradient2(low = 'blue', high = 'red', limits = c(-1,NA),
                       na.value = 'grey80') +
  theme_pubr(legend = 'right') +
  labs(x = '', y = 'Gene')

library(pathview)
## MVA pathway -----------
png.out <- pathview(
  gene.data = v3h.sa.vs.pbs,
  gene.idtype = 'SYMBOL',
  species = 'mouse',
  pathway.id = '00900',
  out.suffix = "mva.path.v3h.sa.vs.pbs",
  kegg.native = T
)

## ROS carcinogenesis pathway ---------
v3h.sa.vs.pbs |>
  pathview(gene.idtype = 'SYMBOL',
           species = 'mouse',
           pathway.id = '05208',
           out.suffix = "ROS.cancer.v3h.sa.vs.pbs")

## apoptosis pathway --------
v3h.sa.vs.pbs |> pathview(
    gene.idtype = 'SYMBOL',
    species = 'mouse',
    pathway.id = '04210',
    out.suffix = "apoptosis.v3h.sa.vs.pbs",
  )

v3l.sa.vs.pbs |> pathview(
  gene.idtype = 'SYMBOL',
  species = 'mouse',
  pathway.id = '04210',
  out.suffix = "apoptosis.v3l.sa.vs.pbs",
)

## apoptosis multispecies pathway -------------
v3h.sa.vs.pbs |> pathview(
  gene.idtype = 'SYMBOL',
  species = 'mouse',
  pathway.id = '04215',
  out.suffix = "apoptosis.mp.v3h.sa.vs.pbs",
  same.layer = F
)

# LISA: UPR TFs -------
savpbs.skin.logfc <-
  read_csv('mission/FPP/zww_sa_mice/results/savpbs.skin.all.logfc.csv')

proj.nm <- 'mission/FPP/zww_sa_mice/'

upr_tf <- c('ATF4','ATF6','XBP1','SREBF1','SREBF2')

upr_target <- read_tsv('mission/FPP/thapsigargin_upr/upr_3tf_targets.tsv')

savpbs_upr <- upr_target |>
  pivot_longer(1:3, names_to = 'TF', values_to = 'gene',
               values_transform = str_to_title) |>
  inner_join(savpbs.skin.logfc, relationship = "many-to-many")

savpbs_upr |>
  ggplot(aes(gene, cluster, fill = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point(shape = 21) +
  facet_wrap(~TF, scales = 'free_x') +
  scale_fill_distiller(palette = 'RdYlBu',
                       values = pretty_distiller(savpbs_upr$avg_log2FC)) +
  theme_pubr() +
  rotate_x_text() +
  labs(title = 'UPR-associated TF downstream genes in SA vs PBS skin')

savpbs.skin.logfc$cluster |>
  unique() |>
  walk(\(x)savpbs.skin.logfc |>
         filter(p_val_adj < .05, avg_log2FC > 1, cluster == x) |>
         select(gene) |>
         write_source_csv(str_c('SAvsPBS_', x), col_names = F))

kera_infct_deg <-
  read_csv('mission/FPP/zww_sa_mice/results/kera.infect.deg.csv')

kera_infct_deg |>
  filter(avg_log2FC > 1, p_val_adj < .05) |>
  select(gene) |>
  write_source_csv('SAvsPBS_Keratinocytes', col_names = F)

c('DC', 'Endothelial cells') |>
  map(\(x)savpbs.skin.logfc |>
         filter(p_val_adj < .1, avg_log2FC > 0, cluster == x) |>
         select(gene) |>
         write_source_csv(str_c('SAvsPBS_', x), col_names = F))

lisa_savpbs <-
list.files('mission/FPP/zww_sa_mice/results/', 'motif_SAvsPBS_',
           full.names = T) |>
  read_tsv(id = 'file') |>
  mutate(file = str_remove_all(file, '.+PBS_|.csv.+'))

lisa_savpbs |>
  filter(factor %in% c('SREBF1','SREBF2')) |>
  mutate(file = fct_relevel(file, 'Keratinocytes', 'V3-hi KC', 'V3-lo KC')) |>
  ggplot(aes(file, -log10(summary_p_value), fill = factor)) +
  geom_col(position = 'dodge2') +
  labs(x = 'Cell type', fill = 'Transcription Factor',
       y = '-log10(p_value)',
       title = 'MVA-associated TF activity upregulated in SA vs PBS skin') +
  geom_hline(yintercept = -log10(.01), linetype = 'dashed') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub +
  rotate_x_text(45)

publish_source_plot('SREBF.SAvsPBS.skin.celltype.LISA', width = 80)
  
lisa_savpbs |>
  filter(str_detect(factor, 'SREBF')) |>
  mutate(file = fct_reorder(file, summary_p_value, .fun = min),
         factor = fct_reorder(factor, summary_p_value, .fun = min)) |>
  ggplot(aes(file, -log10(summary_p_value), fill = factor)) +
  geom_col(position = 'dodge2') +
  labs(x = 'Cell type', fill = 'Transcription Factor',
       y = '-log10(p_value)',
       title = 'UPR-associated TF activity upregulated in SA vs PBS skin') +
  geom_hline(yintercept = -log10(.05), linetype = 'dashed') +
  theme_jpub +
  rotate_x_text(45)

## ctrl V3h vs V3l ---------
ct_v3hvl_deg <-
  read_csv('mission/FPP/zww_sa_mice/results/mice_ctrl_v3h-v3l_deg.csv')

ct_v3hvl_deg |>
  filter(avg_log2FC > 1, p_val_adj < .05) |>
  select(gene) |>
  write_source_csv('PBS_V3hi_vs_V3lo', col_names = F)

pbs_v3hvl_lisa <-
read_tsv('mission/FPP/zww_sa_mice/results/basic_motif_PBS_V3hi_vs_V3lo.csv.lisa.tsv')

pbs_v3hvl_lisa |>
  filter(factor %in% upr_tf) |>
  ggplot(aes('V3h', -log10(summary_p_value), fill = factor)) +
  geom_col(position = 'dodge2') +
  labs(x = 'Cell type', fill = 'Transcription Factor',
       y = '-log10(p_value)',
       title = 'UPR-associated TF activity upregulated in skin cell types') +
  geom_hline(yintercept = -log10(.05), linetype = 'dashed') +
  theme_jpub +
  rotate_x_text(45)

## SA all marker ----------
sa_skin_mrk <-
  read_csv('mission/FPP/zww_sa_mice/results/sa.skin.all.marker.csv')

sa_skin_mrk$cluster |>
  unique() |>
  walk(\(x)sa_skin_mrk |>
         filter(p_val_adj < .05, avg_log2FC > 1, cluster == x) |>
         select(gene) |>
         write_source_csv(str_c('SA_allmarker_', x), col_names = F))

lisa_sa_all <-
  list.files('mission/FPP/zww_sa_mice/results/', 'motif_SA_allmarker',
             full.names = T) |>
  read_tsv(id = 'file') |>
  mutate(file = str_remove_all(file, '.+allmarker_|.csv.+'))

lisa_sa_all |>
  filter(factor %in% upr_tf) |>
  ggplot(aes(file, -log10(summary_p_value), fill = factor)) +
  geom_col(position = 'dodge2') +
  labs(x = 'Cell type', fill = 'Transcription Factor',
       y = '-log10(p_value)',
       title = 'UPR-associated TF activity upregulated in skin cell types') +
  geom_hline(yintercept = -log10(.05), linetype = 'dashed') +
  theme_jpub +
  rotate_x_text(45)

lisa_sa_all |>
  filter(str_detect(factor, 'SREBF\\d$')) |>
  mutate(file = fct_relevel(file, 'Keratinocytes', 'V3-hi KC', 'V3-lo KC')) |>
  ggplot(aes(file, -log10(summary_p_value), fill = factor)) +
  geom_col(position = 'dodge2') +
  labs(x = 'Cell type', fill = 'Transcription Factor',
       y = '-log10(p_value)',
       title = 'SREBF1/2 activity upregulated in SA-infected mouse skin cell types') +
  geom_hline(yintercept = -log10(.05), linetype = 'dashed') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub +
  rotate_x_text(45)

publish_source_plot('SREBF.SA.skin.celltype.LISA', width = 80)

## PBS all marker -------------
pbs_skin_mrk <-
  read_csv('mission/FPP/zww_sa_mice/results/pbs.skin.all.marker.csv')

pbs_skin_mrk$cluster |>
  unique() |>
  walk(\(x)pbs_skin_mrk |>
         filter(p_val_adj < .05, avg_log2FC > 1, cluster == x) |>
         select(gene) |>
         write_source_csv(str_c('PBS_allmarker_', x), col_names = F))

lisa_pbs_all <-
  list.files('mission/FPP/zww_sa_mice/results/', 'motif_PBS_allmarker',
             full.names = T) |>
  read_tsv(id = 'file') |>
  mutate(file = str_remove_all(file, '.+allmarker_|.csv.+'))

lisa_pbs_all |>
  filter(str_detect(factor, 'SREBF\\d$'),
         str_detect(file, 'Kera|KC|Fibro|Gran|Macro')) |>
  mutate(file = fct_relevel(file, 'Keratinocytes', 'V3-hi KC', 'V3-lo KC')) |>
  ggplot(aes(file, -log10(summary_p_value), fill = factor)) +
  geom_col(position = 'dodge2') +
  labs(x = 'Cell type', fill = 'Transcription Factor',
       y = '-log10(p_value)',
       title = 'SREBF1/2 activity upregulated in uninfected mouse skin cell types') +
  geom_hline(yintercept = -log10(.05), linetype = 'dashed') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub +
  rotate_x_text(45)

publish_source_plot('SREBF.PBS.skin.celltype.LISA', width = 80)

# RcisTarget: UPR TFs ------
library(RcisTarget)

data("motifAnnotations_mgi_v9")
motifRankings <-
  importRankings('~/append-ssd/pyscenic/mm10_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather')

celltype_genelist <- sa_skin_mrk |>
  filter(avg_log2FC > 1, p_val_adj < .05) |>
  select(cluster, gene) |>
  pivot_wider(names_from = cluster, values_from = gene, values_fn = list) |>
  as.list() |>
  list_flatten()

sa_all_cistar <- celltype_genelist |>
  cisTarget(motifRankings, motifAnnotations_mgi_v9)

sa_all_cistar |>
  as_tibble() |>
  filter(str_detect(TF_lowConf, 'Srebf') | str_detect(TF_highConf, 'Srebf'))

sa_all_cistar |>
  as_tibble() |>
  distinct(motif, TF_highConf, TF_lowConf) |>
  DT::datatable()
